scripts/Simlation study/start-simulations-sx2.R

In joakimwallmark/PolyOptimalIRT:

library(TestGardener)

# natMath -----------------------------------------------------------------
U         <- scan("data/NatMath.txt","o")
N         <- length(U) # Number of examinees
Umat      <- as.integer(unlist(stringr::str_split(U,"")))
n         <- length(Umat)/N # Number of items
U         <- matrix(Umat,N,n,byrow=TRUE)

#  data in score mode, convert to index mode
U = U + 1

scoreList <- vector("list",32) # option scores
scoreList[[ 1]] = c(0, 1)
scoreList[[ 2]] = c(0, 1)
scoreList[[ 3]] = c(0, 1)
scoreList[[ 4]] = c(0, 1)
scoreList[[ 5]] = c(0, 1)
scoreList[[ 6]] = c(0, 1, 2)
scoreList[[ 7]] = c(0, 1)
scoreList[[ 8]] = c(0, 1)
scoreList[[ 9]] = c(0, 1)
scoreList[[10]] = c(0, 1, 2)
scoreList[[11]] = c(0, 1)
scoreList[[12]] = c(0, 1)
scoreList[[13]] = c(0, 1)
scoreList[[14]] = c(0, 1, 2)
scoreList[[15]] = c(0, 1, 2)
scoreList[[16]] = c(0, 1)
scoreList[[17]] = c(0, 1, 2)
scoreList[[18]] = c(0, 1, 2)
scoreList[[19]] = c(0, 1, 2)
scoreList[[20]] = c(0, 1, 2)
scoreList[[21]] = c(0, 1, 2, 3)
scoreList[[22]] = c(0, 1, 2)
scoreList[[23]] = c(0, 1, 2, 3)
scoreList[[24]] = c(0, 1, 2)
scoreList[[25]] = c(0, 1, 2)
scoreList[[26]] = c(0, 1, 2)
scoreList[[27]] = c(0, 1, 2)
scoreList[[28]] = c(0, 1, 2)
scoreList[[29]] = c(0, 1, 2, 3)
scoreList[[30]] = c(0, 1, 2)
scoreList[[31]] = c(0, 1, 2)
scoreList[[32]] = c(0, 1, 2, 3, 4)
maxScore <- sum(sapply(scoreList, max))
optList <- list(itemLab=NULL, optLab=NULL, optScr=scoreList)

cyc <- c(5, 10)
noBasis <- c(15) # 12 becomes default for natMath data
noBin <- c(25) # 22 becomes default for natMath data
noBasis <- c(7, 9, 12) # 12 becomes default for natMath data
noBin <- c(15, 22) # 22 becomes default for natMath data

## cycle through scenarios
totScenarios <- (length(cyc)*length(noBasis)*length(noBin))
scenario <- 1
set.seed(1233)
start <- proc.time()
for (cy in cyc) {
    for (bin in noBin) {
        for (bas in noBasis) {
            dataList <- make.dataList(U, NULL, optList, scrrng = c(0,maxScore), NumBasis = bas, nbin = bin)
            TGmod <- Analyze(dataList$percntrnk, dataList$thetaQnt, dataList, ncycle=cy, itdisp=FALSE)
            sx2 <- itemFitTestGardener(TGmod, data = U-1, kernelDensity = T)$SX2
            save(sx2, file = paste("data/sx2/natmath ", "kernel",
                                   " cyc", cy,
                                   " basis", bas,
                                   " bins", bin, ".rda", sep = ""))
            sx2 <- itemFitTestGardener(TGmod, data = U-1, kernelDensity = F)$SX2
            save(sx2, file = paste("data/sx2/natmath",
                                   " cyc", cy,
                                   " basis", bas,
                                   " bins", bin, ".rda", sep = ""))
            proc.time()-start
            scenario/totScenarios
            scenario <- scenario+1
        }
    }
}

# hads anxiety --------------------------------------------------------------------
U <- scan("data/hads.txt","o")
N         <- length(U) # Number of examinees
Umat      <- as.integer(unlist(stringr::str_split(U,"")))
n         <- length(Umat)/N # Number of items
U         <- matrix(Umat,N,n,byrow=TRUE)

U <- U[, c(1, 3, 5, 7, 9, 11, 13)]
n <- 7
U <- U + 1

scoreList <- list() # option scores
for (item in 1:7){
    scoreList[[item]] <- c(0:3)
}
maxScore <- sum(sapply(scoreList, max))
optList <- list(itemLab=NULL, optLab=NULL, optScr=scoreList)

cyc <- c(5, 10)
noBasis <- c(10, 12) # 8 becomes default for hads data
noBin <- c(18, 22) # 14 becomes default for 90% of hads data
noBasis <- c(6, 8, 10) # 8 becomes default for hads data
noBin <- c(12, 14) # 14 becomes default for 90% of hads data

## cycle through scenarios
totScenarios <- (length(cyc)*length(noBasis)*length(noBin))
scenario <- 1
set.seed(1233)
start <- proc.time()
for (cy in cyc) {
    for (bin in noBin) {
        for (bas in noBasis) {
            dataList <- make.dataList(U, NULL, optList, scrrng = c(0,maxScore), NumBasis = bas, nbin = bin)
            TGmod <- Analyze(dataList$percntrnk, dataList$thetaQnt, dataList, ncycle=cy, itdisp=FALSE)
            sx2 <- itemFitTestGardener(TGmod, data = U-1, kernelDensity = T)$SX2
            save(sx2, file = paste("data/sx2/hads anxiety kernel",
                                   " cyc", cy,
                                   " basis", bas,
                                   " bins", bin, ".rda", sep = ""))
            sx2 <- itemFitTestGardener(TGmod, data = U-1, kernelDensity = F)$SX2
            save(sx2, file = paste("data/sx2/hads anxiety",
                                   " cyc", cy,
                                   " basis", bas,
                                   " bins", bin, ".rda", sep = ""))
            proc.time()-start
            scenario/totScenarios
            scenario <- scenario+1
        }
    }
}

# hads depression --------------------------------------------------------------------
U <- scan("data/hads.txt","o")
N         <- length(U) # Number of examinees
Umat      <- as.integer(unlist(stringr::str_split(U,"")))
n         <- length(Umat)/N # Number of items
U         <- matrix(Umat,N,n,byrow=TRUE)

U <- U[, c(2, 4, 6, 8, 10, 12, 14)]
n <- 7

#  drop cases containing missing data
keepindex <- U[,7] < 4
sum(keepindex) #  N = 807
U <- U[keepindex,]
U <- U+1

scoreList <- list() # option scores
for (item in 1:7){
    scoreList[[item]] <- c(0:3)
}
maxScore <- sum(sapply(scoreList, max))
optList <- list(itemLab=NULL, optLab=NULL, optScr=scoreList)

cyc <- c(10)
noBasis <- c(10, 12) # 8 becomes default for hads data
noBin <- c(18, 22) # 14 becomes default for 90% of hads data
noBasis <- c(6, 8, 10) # 8 becomes default for hads data
noBin <- c(12, 14) # 14 becomes default for 90% of hads data

## cycle through scenarios
totScenarios <- (length(cyc)*length(noBasis)*length(noBin))
scenario <- 1
set.seed(1233)
start <- proc.time()
for (cy in cyc) {
    for (bin in noBin) {
        for (bas in noBasis) {
            dataList <- make.dataList(U, NULL, optList, scrrng = c(0,maxScore), NumBasis = bas, nbin = bin)
            TGmod <- Analyze(dataList$percntrnk, dataList$thetaQnt, dataList, ncycle=cy, itdisp=FALSE)
            sx2 <- itemFitTestGardener(TGmod, data = U-1, kernelDensity = T)$SX2
            save(sx2, file = paste("data/sx2/hads depression kernel",
                                   " cyc", cy,
                                   " basis", bas,
                                   " bins", bin, ".rda", sep = ""))
            sx2 <- itemFitTestGardener(TGmod, data = U-1, kernelDensity = F)$SX2
            save(sx2, file = paste("data/sx2/hads depression",
                                   " cyc", cy,
                                   " basis", bas,
                                   " bins", bin, ".rda", sep = ""))
            proc.time()-start
            scenario/totScenarios
            scenario <- scenario+1
        }
    }
}